function [SIGMA] = DrawIW(yraw,S_prior,df_prior)

    % S = Scale
    % df = Degrees of Freedom

%   Get initial dimensions of dependent variable
    [T, n] = size(yraw);
  
    % Posterior of SIGMA|ALPHA,Data ~ iW(inv(S_post),df_post)
    df_post = T + df_prior;
    S_post = S_prior + (yraw)'*(yraw);
    SIGMA = inv(wish(inv(S_post),df_post));% Draw SIGMA
    

    
%     draws_prior = nan(4,4,1000);
%     draws_ols = nan(4,4,1000);
%     draws_posterior = nan(4,4,1000);
%     
%     for i = 1:1000
%                 
%         draws_prior(:,:,i) = inv(wish(inv(S_prior),df_prior));
%         draws_ols(:,:,i) = inv(wish(inv((yraw)'*(yraw)),T));
%         draws_posterior(:,:,i) = inv(wish(inv(S_post),df_post));
%         
%     end
%     
%     count=1;
% 
%     for i = 1:4
%         for j = 1:4
%             subplot(4,4,count)
%             hist([squeeze(draws_ols(i,j,:)) squeeze(draws_posterior(i,j,:))],100);
%             count = count+1;
% 
%         end
%     end
    
end

